Prepare Environment
Prepare Data
Disease Map
Point Pattern Analysis
4.1 Density based
4.2 Distance based
Required packages
# install.packages("sf")
# install.packages("tidyverse") #data wrangling
# install.packages("here") #working directory
# install.packages("janitor")
# install.packages("gtsummary")
# install.packages("DT")
# install.packages("stringr")
# install.packages("readxl")
# install.packages("broom")
# install.packages("mapview")
# install.packages("lubridate")
# install.packages("vctrs")
# install.packages("spatialECO") #for NNI
# install.packages("webshot2") #screenshot webpage
Load packages
library(sf)
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ purrr::%||%() masks base::%||%()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(here)
## here() starts at C:/Users/MY PC/OneDrive - Universiti Sains Malaysia/JKNK Spatial R course/spatial_workshop/drhazlienor.github.io
library(janitor)
##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(gtsummary)
library(DT)
library(stringr)
library(readxl)
library(broom)
library(tmap)
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(mapview)
library(lubridate)
library(maptools)
## Loading required package: sp
## Warning: multiple methods tables found for 'elide'
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired during 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
##
## Attaching package: 'maptools'
##
## The following objects are masked from 'package:sp':
##
## elide, sp2Mondrian
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.univar
## spatstat.univar 3.0-0
## Loading required package: spatstat.geom
## spatstat.geom 3.3-2
## Loading required package: spatstat.random
## spatstat.random 3.3-1
## Loading required package: spatstat.explore
## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
##
## spatstat.explore 3.3-1
## Loading required package: spatstat.model
## Loading required package: rpart
## spatstat.model 3.3-1
## Loading required package: spatstat.linnet
## spatstat.linnet 3.2-1
##
## spatstat 3.1-1
## For an introduction to spatstat, type 'beginner'
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(gridExtra)
##
## Attaching package: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combine
library(sparr)
##
##
## Welcome to
## _____ ___ ____ ____ ____
## / ___// _ \/ _ \/ __ \/ __ \
## \__ \/ ___/ __ / ___/ ___/
## ___/ / / / / / / /\ \/ /\ \
## /____/_/ /_/ /_/_/ \__/ \_\ v2.3-15
##
## - type news(package="sparr") for an overview
## - type help("sparr") for documentation
## - type citation("sparr") for how to cite
library(grid)
##
## Attaching package: 'grid'
##
## The following object is masked from 'package:spatstat.geom':
##
## as.mask
library(spatialEco)
##
## Attaching package: 'spatialEco'
##
## The following object is masked from 'package:gridExtra':
##
## combine
##
## The following objects are masked from 'package:spatstat.geom':
##
## is.empty, quadrats, shift
##
## The following object is masked from 'package:spatstat.data':
##
## ants
##
## The following object is masked from 'package:dplyr':
##
## combine
library(plotly)
##
## Attaching package: 'plotly'
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layout
library(gapminder)
library(webshot2)
Types of data: polygon data, population data, point data
Import Data: Load spatial and attribute data (e.g., shapefiles, CSVs) into R using functions like read.csv(), sf::st_read() for vector data, or raster::raster() for raster data.
Clean Data: Handle missing values, correct data formats, and ensure the spatial data is properly projected using sf::st_transform() to standardize coordinate reference systems (CRS).
Explore Data: Use basic summaries (summary(), head()) and visual checks (plot(), View()) to understand data structure.
read polygon data - kelantan map
kel_map <- st_read(here("Map",
"kelantan.shp"))
## Reading layer `kelantan' from data source
## `C:\Users\MY PC\OneDrive - Universiti Sains Malaysia\JKNK Spatial R course\spatial_workshop\drhazlienor.github.io\Map\kelantan.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 66 features and 6 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 371629.6 ymin: 503028.2 xmax: 519479.6 ymax: 690232.8
## Projected CRS: Kertau (RSO) / RSO Malaya (m)
st_geometry() extracts the geometric part (coordinates) of an sf object, allowing you to view or manipulate the spatial features
st_geometry(kel_map)
## Geometry set for 66 features
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 371629.6 ymin: 503028.2 xmax: 519479.6 ymax: 690232.8
## Projected CRS: Kertau (RSO) / RSO Malaya (m)
## First 5 geometries:
## POLYGON ((485501.8 669698.8, 485717.3 669694.6,...
## POLYGON ((487716.5 665649.5, 487615.4 665445.1,...
## POLYGON ((482744.8 660223.4, 482823.6 660137.8,...
## POLYGON ((486936.9 677358.5, 486990.5 677333.7,...
## POLYGON ((490841.7 668783.4, 490906.1 668691, 4...
read population data per subdistrict/mukim per year. You can get the population data from DOSM.
kel_mukim <- read_xlsx(here ("pop_kel.xlsx"))
kel_mukim %>% datatable()
read linelisting cointaining point data in .xlsx format
kel_lepto <- read_xlsx("leptospirosis.xlsx") %>% clean_names()
glimpse(kel_lepto)
## Rows: 1,106
## Columns: 26
## $ diagnosis <chr> "Leptospirosis", "Leptospir…
## $ notifikasi_no <dbl> 2685725, 2728504, 2739963, …
## $ tahun_daftar <dbl> 2016, 2016, 2016, 2016, 201…
## $ epid_daftar <dbl> 6, 9, 11, 12, 18, 18, 19, 1…
## $ age <dbl> 30, 23, 39, 43, 31, 34, 48,…
## $ alamat_semasa_kejadian <chr> "FELDA ARING", "LADANG U&I …
## $ poskod <dbl> 18300, 18300, 18300, 18300,…
## $ latitude_rso <dbl> 478031, 459494, 441802, 488…
## $ longitude_rso <dbl> 548141, 564966, 547551, 547…
## $ latitude_wgs <dbl> 4.944824, 5.034372, 4.89743…
## $ longitude_wgs <dbl> 102.3491, 102.1477, 101.960…
## $ notifikasi_status <chr> "Daftar Kes", "Daftar Kes",…
## $ race <chr> "Foreigner", "Foreigner", "…
## $ kewarganegaraan <chr> "Bukan Warganegara", "Bukan…
## $ gender <chr> "Male", "Male", "Male", "Ma…
## $ nationality <chr> "INDONESIA", "INDONESIA", "…
## $ klasifikasi_kes <chr> "Sporadic", "Sporadic", "Sp…
## $ cara_pengesanan_kes <chr> "Pasif", "Pasif", "Pasif", …
## $ jenis_rawatan <chr> "Wad Perubatan", "Jabatan K…
## $ daerah <chr> "GUA MUSANG", "GUA MUSANG",…
## $ mukim <chr> "CHIKU", "CHIKU", "GALAS", …
## $ lokaliti <chr> NA, "LADANG U & I, CIKU", "…
## $ diagnosis2 <chr> "LEPTOSPIROSIS", "LEPTOSPIR…
## $ sub_diagnosis <lgl> NA, NA, NA, NA, NA, NA, NA,…
## $ negeri <chr> "KELANTAN", "KELANTAN", "KE…
## $ jenis_import_jangkitan_dalam_negara_negeri <chr> NA, NA, NA, NA, NA, NA, NA,…
convert all character to factor (categorical) to allow better understanding of the data
kel_lepto <- kel_lepto %>%
mutate(across(where(is.character), as.factor))
explore the data
summary(kel_lepto)
## diagnosis notifikasi_no tahun_daftar epid_daftar
## Leptospirosis:1106 Min. :2662339 Min. :2016 Min. : 1.00
## 1st Qu.:3005967 1st Qu.:2016 1st Qu.:14.00
## Median :3510632 Median :2018 Median :28.00
## Mean :4701261 Mean :2018 Mean :27.57
## 3rd Qu.:5042478 3rd Qu.:2020 3rd Qu.:42.00
## Max. :9573673 Max. :2022 Max. :53.00
##
## age alamat_semasa_kejadian poskod
## Min. : 0.42 KG MUKA BUKIT : 3 Min. :11510
## 1st Qu.:15.25 KG PANGLIMA BAYU : 3 1st Qu.:16800
## Median :27.00 RPT SG PAS : 3 Median :18000
## Mean :30.95 FELCRA TERATAK BATU: 2 Mean :17543
## 3rd Qu.:45.00 FELDA ARING : 2 3rd Qu.:18300
## Max. :95.00 KAMPUNG BENDANG : 2 Max. :40200
## (Other) :1091 NA's :151
## latitude_rso longitude_rso latitude_wgs longitude_wgs
## Min. : 6 Min. : 63494 Min. :4.611 Min. :101.4
## 1st Qu.:450171 1st Qu.:601283 1st Qu.:5.374 1st Qu.:102.1
## Median :464957 Median :638018 Median :5.764 Median :102.2
## Mean :459729 Mean :627313 Mean :5.651 Mean :102.1
## 3rd Qu.:472980 3rd Qu.:665421 3rd Qu.:6.013 3rd Qu.:102.3
## Max. :502469 Max. :688574 Max. :6.230 Max. :102.5
## NA's :736 NA's :737
## notifikasi_status race kewarganegaraan
## Abai Notifikasi: 5 Aborigines: 27 Bukan Warganegara: 45
## Daftar Kes :1101 Chinese : 9 Warganegara :1061
## Foreigner : 45
## Indian : 5
## Malay :1016
## Others : 4
##
## gender nationality klasifikasi_kes cara_pengesanan_kes
## Female:367 MALAYSIA :1061 Outbreak: 5 Aktif: 6
## Male :739 INDONESIA : 17 Sporadic:1101 Pasif:1100
## MYANMAR : 9
## THAILAND : 7
## BANGLADESH: 6
## INDIA : 3
## (Other) : 3
## jenis_rawatan daerah mukim
## ICU : 58 GUA MUSANG :220 BATU MENGKEBANG: 89
## Jabatan Kecemasan & Kemalangan:385 KUALA KRAI :173 BERTAM : 89
## Jabatan Pesakit Luar :117 MACHANG :143 GALAS : 79
## Lain-lain : 23 PASIR MAS :130 OLAK JERAM : 67
## Wad Kanak-kanak :122 KOTA BHARU :121 ULU SAT : 62
## Wad Perubatan :401 TANAH MERAH: 99 CHIKU : 52
## (Other) :220 (Other) :668
## lokaliti diagnosis2 sub_diagnosis negeri
## LATA REK : 7 LEPTOSPIROSIS:1106 Mode:logical KELANTAN:1106
## PERIA : 7 NA's:1106
## TANAH PUTIH: 7
## RPT SG PAS : 6
## BATU BALAI : 5
## (Other) :1004
## NA's : 70
## jenis_import_jangkitan_dalam_negara_negeri
## KELANTAN: 3
## NA's :1103
##
##
##
##
##
identify observation with missing coordinates. you can update the missing coordinates in the same linelisting and re-run the command
# Identify rows with missing coordinate data
no_coordinates <- kel_lepto %>%
filter(is.na(latitude_wgs) | is.na(longitude_wgs))
# Print rows with no coordinates
print(no_coordinates)
## # A tibble: 0 × 26
## # ℹ 26 variables: diagnosis <fct>, notifikasi_no <dbl>, tahun_daftar <dbl>,
## # epid_daftar <dbl>, age <dbl>, alamat_semasa_kejadian <fct>, poskod <dbl>,
## # latitude_rso <dbl>, longitude_rso <dbl>, latitude_wgs <dbl>,
## # longitude_wgs <dbl>, notifikasi_status <fct>, race <fct>,
## # kewarganegaraan <fct>, gender <fct>, nationality <fct>,
## # klasifikasi_kes <fct>, cara_pengesanan_kes <fct>, jenis_rawatan <fct>,
## # daerah <fct>, mukim <fct>, lokaliti <fct>, diagnosis2 <fct>, …
or, you can proceed with the analysis by removing the observation with no coordinates
kel_lepto2 <- kel_lepto %>% # save to different name to avoid losing the data
filter(!is.na(latitude_wgs),
!is.na(longitude_wgs))
glimpse(kel_lepto2)
## Rows: 1,106
## Columns: 26
## $ diagnosis <fct> Leptospirosis, Leptospirosi…
## $ notifikasi_no <dbl> 2685725, 2728504, 2739963, …
## $ tahun_daftar <dbl> 2016, 2016, 2016, 2016, 201…
## $ epid_daftar <dbl> 6, 9, 11, 12, 18, 18, 19, 1…
## $ age <dbl> 30, 23, 39, 43, 31, 34, 48,…
## $ alamat_semasa_kejadian <fct> "FELDA ARING", "LADANG U&I …
## $ poskod <dbl> 18300, 18300, 18300, 18300,…
## $ latitude_rso <dbl> 478031, 459494, 441802, 488…
## $ longitude_rso <dbl> 548141, 564966, 547551, 547…
## $ latitude_wgs <dbl> 4.944824, 5.034372, 4.89743…
## $ longitude_wgs <dbl> 102.3491, 102.1477, 101.960…
## $ notifikasi_status <fct> Daftar Kes, Daftar Kes, Daf…
## $ race <fct> Foreigner, Foreigner, Forei…
## $ kewarganegaraan <fct> Bukan Warganegara, Bukan Wa…
## $ gender <fct> Male, Male, Male, Male, Mal…
## $ nationality <fct> INDONESIA, INDONESIA, INDON…
## $ klasifikasi_kes <fct> Sporadic, Sporadic, Sporadi…
## $ cara_pengesanan_kes <fct> Pasif, Pasif, Pasif, Pasif,…
## $ jenis_rawatan <fct> Wad Perubatan, Jabatan Kece…
## $ daerah <fct> GUA MUSANG, GUA MUSANG, GUA…
## $ mukim <fct> CHIKU, CHIKU, GALAS, CHIKU,…
## $ lokaliti <fct> NA, "LADANG U & I, CIKU", "…
## $ diagnosis2 <fct> LEPTOSPIROSIS, LEPTOSPIROSI…
## $ sub_diagnosis <lgl> NA, NA, NA, NA, NA, NA, NA,…
## $ negeri <fct> KELANTAN, KELANTAN, KELANTA…
## $ jenis_import_jangkitan_dalam_negara_negeri <fct> NA, NA, NA, NA, NA, NA, NA,…
We want to match the two data based on MUKIM. Ensure the column containing the mukim is named the same for both data.
Use str() to look at the structure of the data.
str(kel_map)
## Classes 'sf' and 'data.frame': 66 obs. of 7 variables:
## $ NEGERI : chr "KELANTAN" "KELANTAN" "KELANTAN" "KELANTAN" ...
## $ DAERAH : chr "BACHOK" "BACHOK" "BACHOK" "BACHOK" ...
## $ MUKIM : chr "BEKLAM" "GUNONG (GUNONG TIMOR)" "MAHLIGAI" "PERUPOK" ...
## $ LELAKI : int 4859 11100 4564 8777 9227 14140 5863 4929 18064 11645 ...
## $ PEREMPUAN : num 4813 10884 4600 8614 8672 ...
## $ JUM_JANTIN: num 9672 21984 9164 17391 17899 ...
## $ geometry :sfc_POLYGON of length 66; first list element: List of 1
## ..$ : num [1:238, 1:2] 485502 485717 485966 486066 486068 ...
## ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA
## ..- attr(*, "names")= chr [1:6] "NEGERI" "DAERAH" "MUKIM" "LELAKI" ...
in ‘kel_map’ polygon data, the column name is ‘MUKIM’
str(kel_mukim)
## tibble [191 × 8] (S3: tbl_df/tbl/data.frame)
## $ subdistrict: chr [1:191] "BEKLAM" "GUNONG (GUNONG TIMOR)" "MAHLIGAI" "PERUPOK" ...
## $ 2016 : num [1:191] 11800 27100 11000 21100 23100 35400 17200 12500 45000 29800 ...
## $ 2017 : num [1:191] 11800 27100 11000 21100 23100 35400 17200 12500 45000 29800 ...
## $ 2018 : num [1:191] 12000 27500 11100 21400 23700 36200 17800 12800 45900 30500 ...
## $ 2019 : num [1:191] 12200 28000 11300 21700 24200 37000 18400 13000 46800 31300 ...
## $ 2020 : num [1:191] 12400 28400 11400 22000 24800 37800 19000 13300 47800 32000 ...
## $ 2021 : num [1:191] 12400 28400 11400 22000 24800 37800 19000 13300 47800 32000 ...
## $ 2022 : num [1:191] 12400 28400 11400 22000 24800 37800 19000 13300 47800 32000 ...
in ‘kel_mukim’ population data, the column name is ‘subdistrict’. We will rename the ‘subdistrict’ variable to ‘MUKIM’.
kel_mukim <- kel_mukim %>%
rename(MUKIM = subdistrict)
str(kel_lepto2)
## tibble [1,106 × 26] (S3: tbl_df/tbl/data.frame)
## $ diagnosis : Factor w/ 1 level "Leptospirosis": 1 1 1 1 1 1 1 1 1 1 ...
## $ notifikasi_no : num [1:1106] 2685725 2728504 2739963 2754040 2811959 ...
## $ tahun_daftar : num [1:1106] 2016 2016 2016 2016 2016 ...
## $ epid_daftar : num [1:1106] 6 9 11 12 18 18 19 17 21 23 ...
## $ age : num [1:1106] 30 23 39 43 31 34 48 55 48 39 ...
## $ alamat_semasa_kejadian : Factor w/ 1073 levels "110, KESEDAR RENOK BARU",..: 108 689 1053 111 676 690 677 220 130 682 ...
## $ poskod : num [1:1106] 18300 18300 18300 18300 18300 ...
## $ latitude_rso : num [1:1106] 478031 459494 441802 488502 496760 ...
## $ longitude_rso : num [1:1106] 548141 564966 547551 547701 539072 ...
## $ latitude_wgs : num [1:1106] 4.94 5.03 4.9 4.95 4.94 ...
## $ longitude_wgs : num [1:1106] 102 102 102 102 102 ...
## $ notifikasi_status : Factor w/ 2 levels "Abai Notifikasi",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ race : Factor w/ 6 levels "Aborigines","Chinese",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ kewarganegaraan : Factor w/ 2 levels "Bukan Warganegara",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ gender : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 1 2 1 ...
## $ nationality : Factor w/ 8 levels "BANGLADESH","CAMBODIA",..: 4 4 4 4 3 3 3 8 4 6 ...
## $ klasifikasi_kes : Factor w/ 2 levels "Outbreak","Sporadic": 2 2 2 2 2 2 2 2 2 2 ...
## $ cara_pengesanan_kes : Factor w/ 2 levels "Aktif","Pasif": 2 2 2 2 2 2 2 2 2 2 ...
## $ jenis_rawatan : Factor w/ 6 levels "ICU","Jabatan Kecemasan & Kemalangan",..: 6 2 6 6 6 6 6 6 6 6 ...
## $ daerah : Factor w/ 10 levels "BACHOK","GUA MUSANG",..: 2 2 2 2 2 2 2 10 2 2 ...
## $ mukim : Factor w/ 65 levels "ALOR PASIR","BADANG",..: 14 14 16 14 14 7 7 48 14 7 ...
## $ lokaliti : Factor w/ 710 levels "AIR BOL","AIR DERAS",..: NA 421 616 170 NA NA NA 667 323 NA ...
## $ diagnosis2 : Factor w/ 1 level "LEPTOSPIROSIS": 1 1 1 1 1 1 1 1 1 1 ...
## $ sub_diagnosis : logi [1:1106] NA NA NA NA NA NA ...
## $ negeri : Factor w/ 1 level "KELANTAN": 1 1 1 1 1 1 1 1 1 1 ...
## $ jenis_import_jangkitan_dalam_negara_negeri: Factor w/ 1 level "KELANTAN": NA NA NA NA NA NA NA NA NA NA ...
in ‘kel_lepto2’ population data, the column name is ‘mukim’. We will rename the ‘mukim’ variable to ‘MUKIM’.
kel_lepto2 <- kel_lepto2 %>%
rename(MUKIM = mukim)
merge population data to polygon
kel <- merge(kel_map,kel_mukim, by.x="MUKIM", by.y="MUKIM", all.x=T, sort=F)
dim(kel)
## [1] 66 14
class(kel)
## [1] "sf" "data.frame"
st_crs(kel)
## Coordinate Reference System:
## User input: Kertau (RSO) / RSO Malaya (m)
## wkt:
## PROJCRS["Kertau (RSO) / RSO Malaya (m)",
## BASEGEOGCRS["Kertau (RSO)",
## DATUM["Kertau (RSO)",
## ELLIPSOID["Everest 1830 (RSO 1969)",6377295.664,300.8017,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4751]],
## CONVERSION["Rectified Skew Orthomorphic Malaya Grid (metres)",
## METHOD["Hotine Oblique Mercator (variant A)",
## ID["EPSG",9812]],
## PARAMETER["Latitude of projection centre",4,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8811]],
## PARAMETER["Longitude of projection centre",102.25,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8812]],
## PARAMETER["Azimuth of initial line",323.0257905,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8813]],
## PARAMETER["Angle from Rectified to Skew Grid",323.130102361111,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8814]],
## PARAMETER["Scale factor on initial line",0.99984,
## SCALEUNIT["unity",1],
## ID["EPSG",8815]],
## PARAMETER["False easting",804670.24,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["Malaysia - West Malaysia onshore."],
## BBOX[1.21,99.59,6.72,104.6]],
## ID["EPSG",3168]]
convert all point data to geometry (sf object). for WGS84, the CRS code is 4326.
lepto_wgs <- st_as_sf(kel_lepto2,
coords = c("longitude_wgs", "latitude_wgs"),
crs = 4326)
confirm CRS is WGS 84
st_crs(lepto_wgs)
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
covert CRS to RSO (because the CRS for polygon map is in RSO-code 3168)
lepto_rso <- st_transform(lepto_wgs, 3168)
plot the map to look for outliers
ggplot() +
geom_sf(data = lepto_rso) +
theme_bw()
if there are outliers, recheck the coordinates, update in the same linelisting and re-run the codes from beginning. or, you can proceed with analysis by excluding points outside the map’s boundaries.
select point only in Kelantan
lepto<- lepto_rso %>%
mutate(within_kel_map = lengths(st_within(lepto_rso, kel_map)))
lepto <- lepto %>%
filter(within_kel_map == 1)
glimpse(lepto)
## Rows: 1,106
## Columns: 26
## $ diagnosis <fct> Leptospirosis, Leptospirosi…
## $ notifikasi_no <dbl> 2685725, 2728504, 2739963, …
## $ tahun_daftar <dbl> 2016, 2016, 2016, 2016, 201…
## $ epid_daftar <dbl> 6, 9, 11, 12, 18, 18, 19, 1…
## $ age <dbl> 30, 23, 39, 43, 31, 34, 48,…
## $ alamat_semasa_kejadian <fct> "FELDA ARING", "LADANG U&I …
## $ poskod <dbl> 18300, 18300, 18300, 18300,…
## $ latitude_rso <dbl> 478031, 459494, 441802, 488…
## $ longitude_rso <dbl> 548141, 564966, 547551, 547…
## $ notifikasi_status <fct> Daftar Kes, Daftar Kes, Daf…
## $ race <fct> Foreigner, Foreigner, Forei…
## $ kewarganegaraan <fct> Bukan Warganegara, Bukan Wa…
## $ gender <fct> Male, Male, Male, Male, Mal…
## $ nationality <fct> INDONESIA, INDONESIA, INDON…
## $ klasifikasi_kes <fct> Sporadic, Sporadic, Sporadi…
## $ cara_pengesanan_kes <fct> Pasif, Pasif, Pasif, Pasif,…
## $ jenis_rawatan <fct> Wad Perubatan, Jabatan Kece…
## $ daerah <fct> GUA MUSANG, GUA MUSANG, GUA…
## $ MUKIM <fct> CHIKU, CHIKU, GALAS, CHIKU,…
## $ lokaliti <fct> NA, "LADANG U & I, CIKU", "…
## $ diagnosis2 <fct> LEPTOSPIROSIS, LEPTOSPIROSI…
## $ sub_diagnosis <lgl> NA, NA, NA, NA, NA, NA, NA,…
## $ negeri <fct> KELANTAN, KELANTAN, KELANTA…
## $ jenis_import_jangkitan_dalam_negara_negeri <fct> NA, NA, NA, NA, NA, NA, NA,…
## $ geometry <POINT [m]> POINT (484205.4 54689…
## $ within_kel_map <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, …
Map all the points on the Kelantan polygon map
ggplot() +
geom_sf(data = kel_map) +
geom_sf(data = lepto) +
theme_bw()
alternatively, you can plot the points on open street view map
library(leaflet)
library(ggmap)
## ℹ Google's Terms of Service: <https://mapsplatform.google.com>
## Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service/>
## OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles/>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
##
## Attaching package: 'ggmap'
##
##
## The following object is masked from 'package:plotly':
##
## wind
leaflet(data = lepto_wgs) %>%
addTiles() %>% # Add default OpenStreetMap map tiles
addCircleMarkers(color = "red", radius = 1, fillOpacity = 0.7)
Map the disease by year.
Adjust the size of the plot using fig.width and fig.height. Adjust the size of the font by increasing or decreasing the text size.
facet_wrap() allows you to stratify the plot based on certain variables (eg: tahun_daftar). Change the variable name if you want to stratify based on other variables (eg: gender)
ggplot() +
geom_sf(data = kel_map) +
geom_sf(data = lepto) +
ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022") +
theme_bw() +
facet_wrap(~ tahun_daftar) +
theme(plot.title = element_text(size = 20), strip.text = element_text(size = 20), axis.text.x=element_text(size=10), axis.text.y=element_text(size=10))
if you want to stratify based on more than one stratified variables (eg: gender and tahun_daftar), use facet_grid()
ggplot() +
geom_sf(data = kel_map) +
geom_sf(data = lepto) +
ggtitle("Map of Leptospirosis Cases in Kelantan by gender for 2016-2022") +
theme_bw() +
facet_grid(gender ~ tahun_daftar) +
theme(plot.title = element_text(size = 24), strip.text = element_text(size = 20), axis.text.x=element_text(size=10), axis.text.y=element_text(size=10))
Interactive plot using plotly
library(plotly)
# Example ggplot object without facet_wrap
p <- ggplot() +
geom_sf(data = kel_map) +
geom_sf(data = lepto) +
ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022") +
theme_bw() +
theme(
plot.title = element_text(size = 20),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10)
)
# Convert the ggplot object to plotly
ggplotly(p)
Interactive map according to year using shiny
library(shiny)
# Define UI for the application
ui <- fluidPage(
titlePanel("Interactive Map of Leptospirosis Cases in Kelantan"),
# Sidebar with two dropdowns: one for year, one for daerah
sidebarLayout(
sidebarPanel(
# Dropdown to select year
selectInput("year", "Select Year:", choices = unique(lepto$tahun_daftar)),
# Dropdown to select daerah
selectInput("daerah", "Select Daerah:", choices = unique(lepto$daerah))
),
# Display the interactive plot
mainPanel(
plotlyOutput("leptoPlot")
)
)
)
# Define server logic
server <- function(input, output) {
output$leptoPlot <- renderPlotly({
# Filter data based on the selected year and daerah
filtered_lepto <- lepto[lepto$tahun_daftar == input$year & lepto$daerah == input$daerah, ]
# Create the ggplot
p <- ggplot() +
geom_sf(data = kel_map) +
geom_sf(data = filtered_lepto) +
ggtitle(paste("Map of Leptospirosis Cases in", input$daerah, "for", input$year)) +
theme_bw() +
theme(
plot.title = element_text(size = 20),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 10)
)
# Convert ggplot to plotly
ggplotly(p)
})
}
# Run the application
shinyApp(ui = ui, server = server)
The points need to be converted to spatial data
As we are goint to analyze the disease for each year, lets extract the leptospirosis cases by year
lepto16 <- subset(lepto, tahun_daftar=="2016")
lepto17 <- subset(lepto, tahun_daftar=="2017")
lepto18 <- subset(lepto, tahun_daftar=="2018")
lepto19 <- subset(lepto, tahun_daftar=="2019")
lepto20 <- subset(lepto, tahun_daftar=="2020")
lepto21 <- subset(lepto, tahun_daftar=="2021")
lepto22 <- subset(lepto, tahun_daftar=="2022")
set observation window
kel_map.owin <- as.owin(kel_map)
plot(kel_map.owin)
Convert to spatial object
lepto16.sp <- as(lepto16, "Spatial")
lepto17.sp <- as(lepto17, "Spatial")
lepto18.sp <- as(lepto18, "Spatial")
lepto19.sp <- as(lepto19, "Spatial")
lepto20.sp <- as(lepto20, "Spatial")
lepto21.sp <- as(lepto21, "Spatial")
lepto22.sp <- as(lepto22, "Spatial")
Convert to planar point pattern (ppp) object
coords16 <- coordinates(lepto16.sp)
lepto16.ppp <- ppp(x = coords16[,1], y = coords16[,2], window = kel_map.owin)
## Warning: data contain duplicated points
coords17 <- coordinates(lepto17.sp)
lepto17.ppp <- ppp(x = coords17[,1], y = coords17[,2], window = kel_map.owin)
coords18 <- coordinates(lepto18.sp)
lepto18.ppp <- ppp(x = coords18[,1], y = coords18[,2], window = kel_map.owin)
## Warning: data contain duplicated points
coords19 <- coordinates(lepto19.sp)
lepto19.ppp <- ppp(x = coords19[,1], y = coords19[,2], window = kel_map.owin)
coords20 <- coordinates(lepto20.sp)
lepto20.ppp <- ppp(x = coords20[,1], y = coords20[,2], window = kel_map.owin)
coords21 <- coordinates(lepto21.sp)
lepto21.ppp <- ppp(x = coords21[,1], y = coords21[,2], window = kel_map.owin)
coords22 <- coordinates(lepto22.sp)
lepto22.ppp <- ppp(x = coords22[,1], y = coords22[,2], window = kel_map.owin)
## Warning: data contain duplicated points
par( mfrow= c(2,4) ) #combine all plot in 1 view, 2 row, 4 column
# 2016
quadr_count_lepto16 <- quadratcount(lepto16.ppp,
nx = 10,
ny = 14)
plot(lepto16.ppp, pch = 20, cex = 0.1, main = "2016")
plot(quadr_count_lepto16, add = TRUE, cex = 2)
# 2017
quadr_count_lepto17 <- quadratcount(lepto17.ppp,
nx = 10,
ny = 14)
plot(lepto17.ppp, pch = 20, cex = 0.1, main = "2017")
plot(quadr_count_lepto17, add = TRUE, cex = 2)
# 2018
quadr_count_lepto18 <- quadratcount(lepto18.ppp,
nx = 10,
ny = 14)
plot(lepto18.ppp, pch = 20, cex = 0.1, main = "2018")
plot(quadr_count_lepto18, add = TRUE, cex = 2)
# 2019
quadr_count_lepto19 <- quadratcount(lepto19.ppp,
nx = 10,
ny = 14)
plot(lepto19.ppp, pch = 20, cex = 0.1, main = "2019")
plot(quadr_count_lepto19, add = TRUE, cex = 2)
# 2020
quadr_count_lepto20 <- quadratcount(lepto20.ppp,
nx = 10,
ny = 14)
plot(lepto20.ppp, pch = 20, cex = 0.1, main = "2020")
plot(quadr_count_lepto20, add = TRUE, cex = 2)
# 2021
quadr_count_lepto21 <- quadratcount(lepto21.ppp,
nx = 10,
ny = 14)
plot(lepto21.ppp, pch = 20, cex = 0.1, main = "2021")
plot(quadr_count_lepto21, add = TRUE, cex = 2)
# 2022
quadr_count_lepto22 <- quadratcount(lepto22.ppp,
nx = 10,
ny = 14)
plot(lepto22.ppp, pch = 20, cex = 0.1, main = "2022")
plot(quadr_count_lepto22, add = TRUE, cex = 2)
Chi-squared goodness-of-fit test
chi_lepto16 <- quadrat.test(lepto16.ppp, nx= 10, ny=14)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
chi_lepto17 <- quadrat.test(lepto17.ppp, nx= 10, ny=14)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
chi_lepto18 <- quadrat.test(lepto18.ppp, nx= 10, ny=14)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
chi_lepto19 <- quadrat.test(lepto19.ppp, nx= 10, ny=14)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
chi_lepto20 <- quadrat.test(lepto20.ppp, nx= 10, ny=14)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
chi_lepto21 <- quadrat.test(lepto21.ppp, nx= 10, ny=14)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
chi_lepto22 <- quadrat.test(lepto22.ppp, nx= 10, ny=14)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
Display result
chi_lepto16.df <- data.frame(Dataset = "lepto16.ppp",
TestStatistic = chi_lepto16$statistic,
PValue = chi_lepto16$p.value)
chi_lepto17.df <- data.frame(Dataset = "lepto17.ppp",
TestStatistic = chi_lepto17$statistic,
PValue = chi_lepto17$p.value)
chi_lepto18.df <- data.frame(Dataset = "lepto18.ppp",
TestStatistic = chi_lepto18$statistic,
PValue = chi_lepto18$p.value)
chi_lepto19.df <- data.frame(Dataset = "lepto19.ppp",
TestStatistic = chi_lepto19$statistic,
PValue = chi_lepto19$p.value)
chi_lepto20.df <- data.frame(Dataset = "lepto20.ppp",
TestStatistic = chi_lepto20$statistic,
PValue = chi_lepto20$p.value)
chi_lepto21.df <- data.frame(Dataset = "lepto21.ppp",
TestStatistic = chi_lepto21$statistic,
PValue = chi_lepto21$p.value)
chi_lepto22.df <- data.frame(Dataset = "lepto22.ppp",
TestStatistic = chi_lepto22$statistic,
PValue = chi_lepto22$p.value)
chi_quadlepto <- bind_rows(
data.frame(Dataset = "2016", chi_lepto16.df),
data.frame(Dataset = "2017", chi_lepto17.df),
data.frame(Dataset = "2018", chi_lepto18.df),
data.frame(Dataset = "2019", chi_lepto19.df),
data.frame(Dataset = "2020", chi_lepto20.df),
data.frame(Dataset = "2021", chi_lepto21.df),
data.frame(Dataset = "2022", chi_lepto22.df))
chi_quadlepto
## Dataset Dataset.1 TestStatistic PValue
## X2...1 2016 lepto16.ppp 1171.8180 6.181358e-181
## X2...2 2017 lepto17.ppp 463.4975 3.361411e-47
## X2...3 2018 lepto18.ppp 479.9587 5.079450e-50
## X2...4 2019 lepto19.ppp 312.4394 6.684179e-23
## X2...5 2020 lepto20.ppp 202.7419 2.378868e-08
## X2...6 2021 lepto21.ppp 173.7245 2.476629e-05
## X2...7 2022 lepto22.ppp 446.4550 2.619867e-44
Monte Carlo test if chi-squared test assumption not met
# run monte carlo test
mc_qlepto16 <- quadrat.test(lepto16.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto17 <- quadrat.test(lepto17.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto18 <- quadrat.test(lepto18.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto19 <- quadrat.test(lepto19.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto20 <- quadrat.test(lepto20.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto21 <- quadrat.test(lepto21.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto22 <- quadrat.test(lepto22.ppp, nx= 10, ny=14, method = "MonteCarlo")
# convert to data frame
mc_qlepto16.df <- data.frame(Dataset = "lepto16.ppp",
TestStatistic = mc_qlepto16$statistic,
PValue = mc_qlepto16$p.value)
mc_qlepto17.df <- data.frame(Dataset = "lepto17.ppp",
TestStatistic = mc_qlepto17$statistic,
PValue = mc_qlepto17$p.value)
mc_qlepto18.df <- data.frame(Dataset = "lepto18.ppp",
TestStatistic = mc_qlepto18$statistic,
PValue = mc_qlepto18$p.value)
mc_qlepto19.df <- data.frame(Dataset = "lepto19.ppp",
TestStatistic = mc_qlepto19$statistic,
PValue = mc_qlepto19$p.value)
mc_qlepto20.df <- data.frame(Dataset = "lepto20.ppp",
TestStatistic = mc_qlepto20$statistic,
PValue = mc_qlepto20$p.value)
mc_qlepto21.df <- data.frame(Dataset = "lepto21.ppp",
TestStatistic = mc_qlepto21$statistic,
PValue = mc_qlepto21$p.value)
mc_qlepto22.df <- data.frame(Dataset = "lepto22.ppp",
TestStatistic = mc_qlepto22$statistic,
PValue = mc_qlepto22$p.value)
# combine rows
mc_qlepto <- bind_rows(
data.frame(Dataset = "2016", mc_qlepto16.df),
data.frame(Dataset = "2017", mc_qlepto17.df),
data.frame(Dataset = "2018", mc_qlepto18.df),
data.frame(Dataset = "2019", mc_qlepto19.df),
data.frame(Dataset = "2020", mc_qlepto20.df),
data.frame(Dataset = "2021", mc_qlepto21.df),
data.frame(Dataset = "2022", mc_qlepto22.df))
mc_qlepto
## Dataset Dataset.1 TestStatistic PValue
## X2...1 2016 lepto16.ppp 1171.8180 0.002
## X2...2 2017 lepto17.ppp 463.4975 0.003
## X2...3 2018 lepto18.ppp 479.9587 0.004
## X2...4 2019 lepto19.ppp 312.4394 0.005
## X2...5 2020 lepto20.ppp 202.7419 0.010
## X2...6 2021 lepto21.ppp 173.7245 0.033
## X2...7 2022 lepto22.ppp 446.4550 0.007
A significant (p-value<0.05) chi-square or Monte-Carlo test result would indicate that the points are clustered and non-randomly distributed, suggesting the presence of spatial processes such as spatial contagion, spatial dependence, or spatial interaction.
par( mfrow= c(2,4) )
inten_lepto16 <-intensity(quadr_count_lepto16)
plot(intensity(quadr_count_lepto16, image = TRUE), main = "2016", las = 1)
plot(lepto16, pch = 20, cex = 0.6, add = TRUE)
## Warning in plot.sf(lepto16, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute
inten_lepto17 <-intensity(quadr_count_lepto17)
plot(intensity(quadr_count_lepto17, image = TRUE), main = "2017", las = 1)
plot(lepto17, pch = 20, cex = 0.6, add = TRUE)
## Warning in plot.sf(lepto17, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute
inten_lepto18 <-intensity(quadr_count_lepto18)
plot(intensity(quadr_count_lepto18, image = TRUE), main = "2018", las = 1)
plot(lepto18, pch = 20, cex = 0.6, add = TRUE)
## Warning in plot.sf(lepto18, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute
inten_lepto19 <-intensity(quadr_count_lepto19)
plot(intensity(quadr_count_lepto19, image = TRUE), main = "2019", las = 1)
plot(lepto19, pch = 20, cex = 0.6, add = TRUE)
## Warning in plot.sf(lepto19, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute
inten_lepto20 <-intensity(quadr_count_lepto20)
plot(intensity(quadr_count_lepto20, image = TRUE), main = "2020", las = 1)
plot(lepto20, pch = 20, cex = 0.6, add = TRUE)
## Warning in plot.sf(lepto20, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute
inten_lepto21 <-intensity(quadr_count_lepto21)
plot(intensity(quadr_count_lepto21, image = TRUE), main = "2021", las = 1)
plot(lepto21, pch = 20, cex = 0.6, add = TRUE)
## Warning in plot.sf(lepto21, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute
inten_lepto22 <-intensity(quadr_count_lepto22)
plot(intensity(quadr_count_lepto22, image = TRUE), main = "2022", las = 1)
plot(lepto22, pch = 20, cex = 0.6, add = TRUE)
## Warning in plot.sf(lepto22, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute
mtext("Intensity Maps of Leptospirosis Cases in Kelantan 2016-2022", side = 1, line = -1, cex = 1.5, outer = TRUE)
###Kernel Density Estimation (KDE)
pre-determined bandwidth (e.g. 5000metres)
kde.lepto16 <- density(lepto16.ppp, sigma = 5000) #2016
kde.lepto17 <- density(lepto17.ppp, sigma = 5000) #2017
kde.lepto18 <- density(lepto18.ppp, sigma = 5000) #2018
kde.lepto19 <- density(lepto19.ppp, sigma = 5000) #2019
kde.lepto20 <- density(lepto20.ppp, sigma = 5000) #2020
kde.lepto21 <- density(lepto21.ppp, sigma = 5000) #2021
kde.lepto22 <- density(lepto22.ppp, sigma = 5000) #2022
par( mfrow= c(2,4) )
plot(kde.lepto16, main = "2016", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto16, add = TRUE)
plot(kde.lepto17, main = "2017", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto17, add = TRUE)
plot(kde.lepto18, main = "2018", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto18, add = TRUE)
plot(kde.lepto19, main = "2019", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto19, add = TRUE)
plot(kde.lepto20, main = "2020", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto20, add = TRUE)
plot(kde.lepto21, main = "2021", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto21, add = TRUE)
plot(kde.lepto22, main = "2022", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto22, add = TRUE)
mtext("Kernel Density Estimate (KDE) Heatmaps of Leptospirosis Cases in Kelantan 2016-2022", side = 1, line = -1, cex = 2, outer = TRUE)
automated bandwidth KDE using bw selector Likelihood Cross Validation
kde.lepto16bw <- density(lepto16.ppp, sigma = bw.ppl(lepto16.ppp))
kde.lepto17bw <- density(lepto17.ppp, sigma = bw.ppl(lepto17.ppp))
kde.lepto18bw <- density(lepto18.ppp, sigma = bw.ppl(lepto18.ppp))
kde.lepto19bw <- density(lepto19.ppp, sigma = bw.ppl(lepto19.ppp))
kde.lepto20bw <- density(lepto20.ppp, sigma = bw.ppl(lepto20.ppp))
kde.lepto21bw <- density(lepto21.ppp, sigma = bw.ppl(lepto21.ppp))
kde.lepto22bw <- density(lepto22.ppp, sigma = bw.ppl(lepto22.ppp))
par( mfrow= c(2,4) )
plot(kde.lepto16bw, main = "2016", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto16bw, add = TRUE)
plot(kde.lepto17bw, main = "2017", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto17bw, add = TRUE)
plot(kde.lepto18bw, main = "2018", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto18bw, add = TRUE)
plot(kde.lepto19bw, main = "2019", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto19bw, add = TRUE)
plot(kde.lepto20bw, main = "2020", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto20bw, add = TRUE)
plot(kde.lepto21bw, main = "2021", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto21bw, add = TRUE)
plot(kde.lepto22bw, main = "2022", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto22bw, add = TRUE)
mtext("Kernel Density Estimate (KDE) Heatmaps of Leptospirosis Cases in Kelantan 2016-2022", side = 1, line = -1, cex = 2, outer = TRUE)
Peaks: Peaks in the KDE plot represent areas where data points are concentrated, indicating higher density or probability of observations in those regions. The height of a peak shows how likely it is to find data points within that range.
Spread/Width: The width or spread of the KDE curve indicates the variability of the data. A wider curve suggests that the data points are more spread out, while a narrower curve indicates that the data points are concentrated in a smaller range.
Shape: The shape of the KDE plot can give insights into the distribution of your data. For instance:
Unimodal: A single peak suggests the data follows a normal-like distribution or another unimodal distribution.
Bimodal/Multimodal: Multiple peaks indicate that the data may have more than one mode, or there are distinct groups or clusters in the dataset.
The lighter the color, the higher the density
The Average Nearest Neighbor (ANN) method measures the distance from each point in a point pattern to its nearest neighboring point and calculates the average of these distances.
This method focuses on local point-to-point relationships and is used to detect clustering or dispersion at a single scale (i.e., based on the nearest neighbor distance).
ann16 <- round(mean(nndist(lepto16.ppp)), 2)
ann17 <- round(mean(nndist(lepto17.ppp)), 2)
ann18 <- round(mean(nndist(lepto18.ppp)), 2)
ann19 <- round(mean(nndist(lepto19.ppp)), 2)
ann20 <- round(mean(nndist(lepto20.ppp)), 2)
ann21 <- round(mean(nndist(lepto21.ppp)), 2)
ann22 <- round(mean(nndist(lepto22.ppp)), 2)
annual_ann <- data.frame(
Year = c(2016, 2017, 2018, 2019, 2020, 2021, 2022),
Nearest_Neighbor_Distance = c(ann16, ann17, ann18, ann19, ann20, ann21, ann22)
)
print(annual_ann)
## Year Nearest_Neighbor_Distance
## 1 2016 1793.08
## 2 2017 3282.28
## 3 2018 2792.46
## 4 2019 4385.30
## 5 2020 7358.82
## 6 2021 6394.20
## 7 2022 2533.81
If the index (average nearest neighbor ratio) is less than 1, the pattern exhibits clustering. If the index is greater than 1, the trend is toward dispersion.
nni16 <-nni(lepto16)
## Warning: data contain duplicated points
nni17 <-nni(lepto17)
nni18 <-nni(lepto18)
## Warning: data contain duplicated points
nni19 <-nni(lepto19)
nni20 <-nni(lepto20)
nni21 <-nni(lepto21)
nni22 <-nni(lepto22)
## Warning: data contain duplicated points
annual_nni <- bind_rows(
data.frame(Dataset = "2016", nni16),
data.frame(Dataset = "2017", nni17),
data.frame(Dataset = "2018", nni18),
data.frame(Dataset = "2019", nni19),
data.frame(Dataset = "2020", nni20),
data.frame(Dataset = "2021", nni21),
data.frame(Dataset = "2022", nni22)
)
annual_nni
## Dataset NNI z.score p expected.mean.distance
## 1 2016 0.5974538 -14.6318827 1.758528e-48 3001.196
## 2 2017 0.8069370 -4.6425708 3.441007e-06 4067.577
## 3 2018 0.6811086 -8.1163418 4.804463e-16 4099.872
## 4 2019 0.7605712 -4.5112102 6.445880e-06 5765.793
## 5 2020 1.0425996 0.5281551 5.973917e-01 7058.142
## 6 2021 0.9495579 -0.5948619 5.519358e-01 6733.872
## 7 2022 0.6918232 -8.9993010 2.271594e-19 3662.505
## observed.mean.distance
## 1 1793.076
## 2 3282.278
## 3 2792.458
## 4 4385.297
## 5 7358.816
## 6 6394.201
## 7 2533.806
Empirical values greater than theoretical (Poisson) values suggest clustering. Envelope denote simulations to test for CSR
par( mfrow= c(2,4) )
G_lepto16 <- plot(envelope(lepto16.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto17 <- plot(envelope(lepto17.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto18 <- plot(envelope(lepto18.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto19 <- plot(envelope(lepto19.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto20 <- plot(envelope(lepto20.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto21 <- plot(envelope(lepto21.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto22 <- plot(envelope(lepto22.ppp, Gest, nsim = 99, verbose = FALSE))
Empirical values greater than theoretical (Poisson) values suggest clustering. Envelope denote simulations to test for CSR Takes long time, for demonstration purpose, we will set the Monte Carlo simulation to only 3
par( mfrow= c(2,4) )
K_lepto16 <- plot(envelope(lepto16.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto17 <- plot(envelope(lepto17.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto18 <- plot(envelope(lepto18.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto19 <- plot(envelope(lepto19.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto20 <- plot(envelope(lepto20.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto21 <- plot(envelope(lepto21.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto22 <- plot(envelope(lepto22.ppp, Kest, nsim = 3, verbose = FALSE))
other summary function - Lest, Fest, pcf